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ABSTRACT . A review is presented of recent developments in the calcula- 
tion or design parameters for fracture mechanics by the method of lines 
(MOL). Three-dimensional elastic and elasto-plastic formulations are exam- 
ined and results are reported from previous and current research activi- 
ties. The application of MOL to the appropriate partial oifferential equa- 
tions of equilibrium leads to coupled sets of simultaneous oroinary differ- 
ential equations. Solutions of these equations are obtained by the Peano- 
Baker and by the recurrance relations methods. The advantages and limita- 
tions of both solution methods from the computational standpoint are sum- 
marized. 

INTRODUCTION . The main goal of fracture mechanics is the prediction of 
the load at which a structure weakened by a crack will fail. Knowledge of 
the stress and displacement fields near the crack tip is of fundamental im- 
portance in evaluating this load at failure. Because of the geometric sin- 
gularity associated with any crack type problem, there is almost no possi- 
bility of a simple closed form type of solution. For this reason, three- 
dimensional elastic solutions have been obtained only for a restricted class 
of problems. Furthermore, the calculation of stress and strain distribu- 
tions in elasto-plastic/work hardening materials containing inherent crack- 
like flaws is a non-linear and three-dimensional problem. Due to the finite 
boundary effect and the nonlinearity of the material response, solutions in 
existence are obtained almost exclusively through numerical computer methods 
of continuum mechanics. Notable among these are the finite element method 
[1,2], the finite difference method l 3], and the boundary integral equation 
method [4]. These methods are useful in solving either elastic or elasto- 
plastic fracture mechanics problems; it is known, however that practical 
problems usually require a very large amount of data storage and computation 
time. 


An alternate semi-analytical method suitable for the solution of crack 
problems is the line method of analysis. Successful application of this 
method to finite geometry solids containing cracks has been demonstratea 
recently for both elastic [5] and elasto-plastic [6] problems. Although the 
concept of the line method for solving partial differential equations is not 
new [7,8j, its application in structural analysis has been limited to simple 
examples [9]. By far the most common approach to fracture problems has been 
the finite element method, and it is the purpose of this paper to review a 
simple, systematic, alternate method, the method of lines (Mul) for these 
problems. 

The line method lies midway between completely analytical and dis- 
cretized numerical methods. The basis of this technique is the substitution 
of finite differences for the derivatives with respect to all the independ- 
ent variables except one for which the derivatives are retained. This ap- 
proach replaces a given partial oifferential equation with a system of ordi- 
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nary differencial equations whose solutions can then be obtained, at least 
in some cases, by analytic methods. These equations describe the dependent 
variable along lines which are parallel to the coordinate in whose direction 
the derivatives were retained. Application of the line method is most use- 
ful when the resulting ordinary differential equations are linear and have 
constant coefficients L 103 . 

Ai .nherent advantage of the line method over other numerical metnoas 
is that good results are obtained from the use of relatively coarse grids. 
This use of a coarse grid is permissible because parts of the solutions are 
obtained in terms of continuous functions. It is known that MOL methods 
tend to keep the advantages ana discard the disadvantages of both the ana- 
lytical and grid methods, thereby leading to accurate solutions with minimun 
computation times. The disadvantage of MOL, on the other hana, is that ic 
tenas to become numerically unstable as the number of dividing lines in- 
creases and the finite difference strip sue becomes too small L&»ll,12j. 

To realize a very fine space discretization with this method would require 
word length with much larger number of bits, leading to excessive require- 
ments on computer resources. Current research emphasis in MOL solution 
methods is to overcome this problem in engineering applications Ll^J* 

G OVERNING EQUATIONS AND MOL FORMULATION . It is assumed, for simplicity 
of tliis presentation, that the material is homogeneous, isotropic and that 
the oeformations are quasi-static and small. The structure is assumed to be 
elastic f irst and the elastic solution is taken to be known before the in- 
cipient loading is applied. As loading gradually increases, the structure 
becomes elasto-plastic and the governing equations are written in terms of 
displacement increments. Using the standard summation convention, thr 
Navier equations for the elastic problem in terms of displacement*. u 1# are 

u i,jj + (rh^) u j,.p (1 

end for the elasto-plastic regime, the displacement increments, ou-j, can 
be obtained from 



where the body forces are assumed to be zero, dc p is the effective plas- 
tic strain increment, is the stress deviator tensor ana a e is the 
equivalent stress. In tne plastic region the von Mises yield condition ana 
the associated Prandtl-keuss flow rule is taken to prevail. The incremental 
stress-strain relations are obtained as l&J. 


- ae U * (r^77) °*kk 6 1 j *l(^)[3(l - 2.)] c kk‘ij 

where v, G, E are the conventional elastic properties, o-jj is tne 
Kronecker delta and oij are the stresses. 





In o.'oer to solve equations (1) or (2), we apply MOL and reduce these 
equations to systems of simultaneous ordinary differential equations. For 
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problems in Cartesian coordinates, the reyion is discretized by x, y and 
z-directional lines as shown in figure 1. The displacements along the 
x-directional lines are defined as uj, u^, . . u r The deriva- 
tives of the y-directional displacements on these lines with respect to y 
are defined as v'li» v 1 j 2 * • • •* v * j £ , anu the derivatives of the 
z-directional displacements with respect to z are defined as w'h, 
w ' J 2 , . . . , w * | £ . When these definitions are used the ordinary differ- 
ential equation along a generic line ij (a double subscript is used here 
for simplicity of writing and the subscripts ooviously are not relateo to 
those in the equations) in figure 1, using central differences with trunca- 
tion errors of 0 (h 2 ), may be written as 



Similar differential equations are obtained along the other x-directiona 1 
lines. The set of ! second order differential equations represented by 
(4) can be reduced to a set of Zi first order differential equations by 
treating the derivatives of the u's as an additional set of e unknowns. 
The resulting equations can now be written as a single first order matrix 
differential equation 


^ - A U + R (x ) (t>> 

where U and R are column matrices of Za elements each and Aj is 
Zi x 2« matrix of coefficients. In a similar manner to solve the other two 
havier equations for the elastic problem, we construct ordinary differential 
equations along the y- and z-directional lines, respectively. These equations 
are also expressed in an analogous form to equation (6); they are 

37-V* sl>) m 
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dW 

dz 


A^W + T(z) 
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Equations (6) to (8) are linear, first-order, ordinary differential equa- 
tions. They are, however, not independent, but are coupled through the vec- 
tors R, S and T. 

Noting that a second order ordinary differential equation can satisfy 
only a total of two boundary conditions and since three-dimensional elastic- 
ity problems have three conditions at every point of the bounding surface, 
the shear stress boundary data must be incorporated into the differential 
equations of the surface lines. The application of the specified shear 
conditions permits the use of a single layer of oounoary image lines when 
surface line differential equations are constructed. 

For an elasto-plastic solid the governing differential equations for 
displacement increments and the incremental stress-displacement relations 
are found in |,6]. The x-directional displacement increments, in an analo- 
gous manner to equation (6), can be obtained from 

(dll) - AjdU + dR(x) 

where the coupling vector dR(x) contains mixed derivative terms for elastic 
and plastic regions in aodition to terms involving the ratio of dep/o e . 

The system of ordinary differential equation (6) can be solved by any 
of a number of standard techniques. The method employed in [5,b,9j is the 
Peano-Baker method of integration. The solution can oe written as 

A.x A,x f x -A,n 

U(x) = e U(0) + e / e R(n) dn UO) 

J 0 

where U(0) is the initial value vector determined from the boundary condi- 
tions and the matrizant e A l x is generally evaluated by its matrix 
series. For larger values of x, when corvergence becomes slow, adoitive 
formulas may be used. In addition, similarity transformations can be used 
to diagonalize the coefficient matrix Ai. It shoulo De noted that, in 
general, the matrix Aj is a function of Poisson's ratio ano tne coordi- 
nate finite difference increments. Uniform line spacing in the three coor- 
dinate directions makes closed form diagonal ization of Aj possible. 

However, refinement of the mesh wit! uniform line spacing rapidly increases 
the required computer time and storage as well as raises the probability of 
numerical difficulties in the matrix exponential power series computations. 
Consequently, variable mesh spacing is recommended as one method of ootain- 
ing improved answers without an increase in problem size. 

Most of the recently obtained MOL solutions in fracture mechanics in- 
volved the use of finite difference formulas with truncation errors of 
0(h2). Current work by Mendelson ano Alam l 13 j, uses higher order finite 
difference approximations as an alternate method of obtaining more accur^., 
results. These five point finite difference approximations for tne first 
and second y-aerivatives of a function f(x,y,z) at (x,y,z) can written 
as [8], 
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Since the evaluation of the matrix exponential pouter series becomes increas- 
ingly difficult with the coefficients obtained through the use of equations 
(11), a recurrence relations method is used to solve the resulting systems 
of ordinary differential equations. An error analysis in L#] indicates that 
approximately six times as many dividing lines must be used with 0(h2) 
approximations to get equivalent accuracy to that obtained when equations 
(11) are used. 


The solution of equation (6) by recurrance relations can be obtained by 
taking N equal intervals along the x-axis, each having a length of h m . 

Then by using finite differences for (du/dx) m and average values of 
(Aju + R) m , the following linear recurrance formula will be obtained 
[13], expressing U m in terms of u m _i : 


where L m is a known function of h m and Aj while will de- 
pend on R in addition to h m and Aj. We can also express in 
terms of U} by repeated application of equation (12), leading to 


(U) 


where tfo and F m are known functions of L m and M m . By suit- 
able partitioning of the D and F matrices at the boundaries, we can use 
given boundary data at the last station, u n , to calculate unknown elements 
of Uj, where n « N + 1. The advantage of using equation (13) to calcu- 
late dm, m » 1, 2, . . . , n, is that the coefficient matrix A^ has no 
limitation on its format or on the arrangement of its elements. The inter- 
val h m can be decreased to any fraction of h x , the initially estab- 
lished finite difference increment obtained from the application of MOL. 
Results of test problems indicate that two or three subintervals are ade- 
quate for the solution of a typical problem. 


All of the MOL work in three-dimensional fracture mechanics has been 
done using double precision arithmetic. With larger word sues, 128 bits or 
greater, improved results can be obtained or more lines can be used. 

Typical problems on third generation computers can usually be handled witn 
100 to 200 lines in each direction using Cartesian coordinates, and up to 20 
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lines in each direction using cylindrical coordinates. Corresponding CPU 
times are of the order of 3 minutes for the elastic solution and 25 minutes 
for the entire elasto-plastic problem (6]. Most of the computer time is 
spent in decoupling the three systems of simultaneous ordinary difrerential 
equations which arise in a general problem. Decoupling is done by a succes- 
sive approximation procedure, with cyclic resubstitution of the obtained 
solutions into tne coupling vectors and the boundary conditions good numer- 
ical convergence behavior was observed. 

Elastic solutions have been obtained for typical fracture test specimen 
geometries such as the central crack, single edge crack, double edge crack 
and rectangular surface crack problems, all with uniform tensile loadings 
normal to the crack plane. In cylindrical coordinates, the problem of an 
embedded penny-shaped crack and the externally cracked circular cylinder 
tension have been treated. Although shear and torsional ly loaded specimens 
have not been analyzed previously by MOL, the necessary boundary conditions 
for these problems can be imposed without difficulties. Presently, the 
thumbnail -crack problem is being investigated in connection with applying 
MOL to curved crack boundaries. In addition, it seems that the common 
three-point bend specimen could also be analyzed in a systematic manner 
using this method. 

Elasto-plastic solutions of certain crack problems have also been ob- 
tained using the incremental displacement formulation in connection with 
MOL. The nonlinear response of finite length cylinaers with external annu- 
lar cracks and a finite thickness rectangular plate containing a through- 
thickness central crack under uniaxial tension were studied in detail. In 
addition to the stresses and displacements, fracture mechanics parameters 
such as the stress intensity factor, the J-integral and the load versus load 
point displacement plots were determined. 

STRESS INTENSITY AND STRAIN ENERGY DENSITY FACTORS . Since the applica- 
tion of boundary conoitions for MOL allows the crack tip to remain between 
two successive node points, the exact location of crack tip together with 
the determination of K values for the elastic case is done by the first 
two terms in the Williams eigenfunction expansion. The crack face displace- 
ment and maximum normal stress near tne crack front are used to find the 
coefficients in the two-term expansion. Assuming that y * U is the crack 
plane and that the crack is under normal tensile loading, we have 


v * oKj 




U4) 
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where o is a function of Poisson's ratio, the stress singularity is 
assumed to be -1/2, r is the crack edge position correction, v 1 tne 
crack opening displacement, oy is the maximum normal stress ano 
and Lj are the mode I Williams expansion coefficients, using c.splacc- 
ment data from three adjacent nodes to the crack edge in equation (14), 
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ues of oK}, Lj/Ki and r are calculated for each value of z, with 
ft also measured from the halfway point between nodes specifying boundary 
stresses and displacements, respectively. Substituting values of 
l i/Kj and r into equation (15), we can calculate K} as a function 
ot the corrected crack edge distance, p ■ R - r. Note that a would be 
equal to 3.56 for the plane strain case and 4.0 for the plane stress case 
with v - 1/3. Another approach to calculate Kj is to first determine 
the 0} Integral and then use the linear elastic - K} relation 
of the form Lb] 



Linear fracture mechanics technology assumes then that the crack will 
propagate If K| reaches Its critical value Kj c , usually called Its 
fracture toughness. It should be noted that the K-concept Is restricted to 
symmetric systems with the applied loads perpendicular to the crack plane 
and the crack propagating In a self-similar manner. In general, a complete 
description of the crack border stress field requires three stress Intensity 
factors, and a mixed mode fracture criterion Is needed to predict failure. 

To this end Sih [14] has defined a strain energy density factor, S, as 

S ■ a n lc i + * a i 2 * 1*2 + a 22 k 2 + *33^3 

where the coefficients a^< (i, j ■ 1, 2) depend on the material constants 
and the angles e and u. Consistent with equation (17), the local 
stresses near the crack tip are of the form £ 14 j, 

o x cos | (l - sin | sin ^ 

y/ip cos u 

1 S in | (2 * cos j cos + 0(1) (lb; 

V2p cos u 

Note that ki « K\ly/lt and the coefficients a^j are then given by 
equations of the form 

166 cos u a^ - (3 - 4v - cos e)(l + cos ©) (19) 

As can be seen from equation (17), the stress intensity factors k-j still 
play an important role in the fracture process. Hence, the correct deter- 
mination of these factors is a necessary step in the safe design of any 
structure. In the strain energy oensity failure criterion, it is assumeo 
that the minimum value of S yields the direction of crack initiation ana 
that the critical value of S^p, S c , determines incipient fracture and 
is an intrinsic material property independent of the loading conditions and 
crack configurations. 

Other multi-mode failure criteria have been proposed previously in the 
literature and a brief description of each fracture theory along with 
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limited experimental data can t>e found in |.i.6,lbj. In general, little dif- 
ference exists between theories predicting mode 1, mode 11 interaction ana 
combined damage cracx propay at ion direction. 

NUMERICAL RESULTS . A great amount of experimental work has been oone 
in fracture mechanics using cracked specimens. Selected results for some 
common specimen geometries are shown in figures 2 to 12. Figure 2 shows a 
rectangular oar unoer normal tensile loading containing a traction free, 
through-thickness, central crack. The crack opening displacement at tne 
middle and at the surface of the bar is plotted in figure 3, along with 
Raju's finite element results. Stress intensity factor variations as a 
function of oar thickness are shown in figure 4. Variation of the con- 
straint parameter a, defined as the quantity o 2 /Lv(ox + oyjj, along 
tne plate thickness is shown in figure 6. Note that for plane strain e - 1 
anc for plane stress case it vanishes. Figure b shows a oar with uniform 
tension containing a rectangular surface crack. Surface crack opening 
displacement as a function of crack depth is shown in figure 7 while the 
variation of the maximum normal stress oy is shown in figure b for a 
selected crack geometry. For the same rectangular surface crack proDiem, a 
plot of Kj along the crack periphery is shown in figure 9. The ois- 
cretuation of an externally cracked cylindrical fracture specimen is shown 
in figure 10. Crack face displacements for various crack lengths are 
plotted in figure 11 while the variation of Kj with cracx length is 
shown in figure 12. It is oovious from these results that a variety of 
plots familiar to the fracture mechanics community can oe constructed, since 
MOL methods give complete fielo solutions. 

CONCLUDING HEhhKKS . Tne line method of analysis is a practical ap- 
pro a cTT7orTnFsinutTon of three-dimensional crack problems, at least for 
bodies with reasonably regular boundaries, both elastic anu inelastic solu- 
tions can oe ootainea. Just how efficient the methoo is or can oe maoe is 
not fully established. It is known, however, that good results are obtaineo 
from the use of relatively coarse grios. Interestingly, emplacements ana 
normal stresses are determined with equal accuracy since numerical differ- 
entiation of the displacements is not required. Applications to curveu 
boundaries, bending or shear modes of loading ana variaole mesh spacing are 
some of the current areas that need additional investigations. Furthermore, 
it seems that MOL could also be useo to stuoy the stable cracx growth be- 
havior of engineering materials. 
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l*> THREE SETS OF LINES PARALLEL TO AND 7- 
COORDINATES AND PERPENDICULAR TO CORRESPONDING 
COORDINATE PUNES 



ft) SET OF INTERIOR LINES PARALLEL TO (-COORDINATE. 
Figure 1 - Sets al lines pertllei to Certesien coordinites 



ill RECTANGULAR BAR WITH THROUGH-THICKNESS CENTRAL 
CRACK. 
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ltd 0ISCRETI2ED REGION Of RECTANGULAR BAR SVITH 
THROUGH-THICKNESS CENTRAL CRACK 
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Figure 5. * Variation of p • o z /[v(o x + o y )] on crack 
plane, y * 0, across plate thickness, (cm) * 0.5174. 
For ideal plane strain case: p * 1, and for plane stress: 
p-a 


RECTANGULAR 
SURFACE CRACK- 






Figure ft H»r with rectangular surface crack under uniform tension 



la) DIMENSIONLESS y-DIRECTIONAL NORMA!. STRESS 
VARIATION ALONG BAP WIDTH. 
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lb) DIMENSIONLESS y-DIRECTlONAL NORMAL STRESS 
VARIATION ACROSS BAR THICKNESS 

Figure 8. - Dimension less y-directionai nori.ia! stress 
distribution in the crack plane for a bar under uniform 
tension containing a rectangular surface crack. 
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